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Abstract. We describe a further development of the stochastic state selection 
method, a new Monte Carlo method we have proposed recently to make numerical 
calculations in large quantum spin systems. Making recursive use of the stochastic state 
selection technique in the Lanczos approach, we estimate the ground state energy of 
the spin-1/2 quantum Heisenberg antiferromagnet on a 48-site triangular lattice. Our 
result for the upper bound of the ground state energy is —0.1833 ± 0.0003 per bond. 
This value, being compatible with values from other work, indicates that our method is 
efficient in calculating energy eigenvalues of frustrated quantum spin systems on large 
lattices. 

PACS numbers: 05.10.Ln,02.70.Ss,75.10.Jm 
1. Introduction 

In numerical studies of quantum spin systems, one of widely used approaches is the 
quantum Monte Carlo method. This method has helped us greatly to understand many 
properties of non-frustrated quantum spin systems, especially of the spin-1/2 quantum 
Heisenberg antiferromagnet on bipartite lattices 12 El H] • Nevertheless, the method is 
ineffective for frustrated systems owing to the so-called sign problem. It is quite difficult, 
therefore, to draw any definite conclusion from numerical calculations of two-dimensional 
large systems of fermions or frustrated spins. Yet there are active studies of numerical 
methods to investigate these systems. One of them is the path-integral renormalization 
group method for fermion systems developed by Kashima and Imada [5 a . Using the 
Slater determinant as the basis state they improve the exact diagonalization method 
with the truncation of the Hilbert space. Another method is proposed by Sorella j^j, 
who extends the fixed node method making full use of insights into the physics of the 
target system. One should also note the work of Henelius and Maeshima et al IE] of 
extending the density matrix renormalization group method introduced by White [9|linj. 

Recently we have developed another Monte Carlo method, which we call the 
stochastic state selection (SSS) method, to calculate eigenvalues in large quantum 
systems [TT] IT2*| EH Ej- The SSS method has little in common with the ordinary 
Monte Carlo methods because it is not based on importance samplings. A new type 
of stochastic algorithm in this method enables us to select a relatively small number of 
elements from a vast vector space in a mathematically justified manner. Using those 
selected elements we calculate inner products. It is guaranteed that we can obtain 
correct values of these inner products through the statistically averaging process. 

So far we have used the SSS method in combination with the power method, since 
it is a simple and straightforward way of applying the SSS method to the numerical 
study on energy eigenvalues. In this paper we employ the Lanczos approach instead of 
the power method. The Lanczos method gives us reliable results on small lattices for 
which we can keep every state in the vector space; for larger lattices we need to make 
some truncations because the vector space becomes huge. Use of truncated states in 
the Lanczos method, however, has been unsuccessful in either theoretically justifying its 
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methodology or in obtaining better numerical eigenvalues ^5] • Our purpose in this work 
is to show that for a large frustrated system we can evaluate the coefficients necessary 
in the Lanczos approach by means of the SSS method. 

As a concrete example, we calculate an energy eigenvalue on a iV s -site triangular 
lattice for the ground state energy of the spin- 1/2 quantum Heisenberg antiferromagnet. 
The Hamiltonian of the system is 



where <Xj denotes Pauli spin matrices on the i-ih site of a triangular lattice with N s 
sites and the sum runs over all iV&(= 3N S ) bonds of the lattice. The coupling J 
is set to 1 throughout this paper. The reason why we study this system here is 
that, as is well-known, it is one of the systems that is strongly frustrated in two 
dimensions |H UHl HZl HH1 HHI 1201 EH 1221 1231 As far as we know, reported results 
on the model by means of the exact diagonalization method are only for lattices smaller 
than or equal to 36 sites JH]. The largest lattice size is also 36 in our previous work |13j . 
where the recursive SSS method and the power method are employed. In the present 
study with an improved approach, we first calculate the lowest energy eigenvalue of the 
N s = 27, S z = 1/2 system in order to confirm that our new approach works well. For 
this system we mimic conditions of a virtual small computer so that we can conclude the 
results obtained in reduced vector space are reliable (table HJ). We then proceed to the 
N s = 48, S z = system (table EJ, for which our result is E/JN b = -0.1833 ± 0.0003. 
A happy combination of the recursive SSS method with the Lanczos approach in this 
paper enables us to estimate the ground state energy of the 48-site system. 

The plan of this paper is as follows. In the next section we explain our method. 
Section El is devoted to a detailed account of the trial state on triangular lattices. 
Calculations with stochastic state selections are described in section HJ The final section 
gives a summary and discussions. 

2. Method 

In this section we present brief descriptions of the SSS method and our Lanczos 



2.1. SSS method 

The stochastic state selection is realized by a number of random variables. Let us 
expand a state | (f>) by some basis {| i)}, | 0) = Y^/iLi Q | i). Then we generate a random 
variable rji following to the on- off probability function 




(1) 



approach. 
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A positive parameter e common to all Pi(i]) (i = 1,2, ■ ■ ■ , N) controls the reduction 
rate. Note that rji = ai or rji = and statistical averages are ((rji)) = 1 and = ai. 
A random choice operator M{ v } is defined by 

N 

M {v} = Y,\j)^{j\ • (3) 

J'=l 

Using this Mr^} we obtain a state M{ v } \ 4>) — S c iVi I > which has fewer non-zero 
elements than | 0). An expectation value ((f> \ O \ <f>) with an operator O is exactly 
equal to the statistical average (( (0 | OM^ \ <fi) )). When we want to emphasize that 
different random choice operators are used, we will denote them by M^^y (k — 1, 2, • • •). 

2.2. Lanczos approach 

In the Lanczos approach we start from a state | ip^) and calculate, for some L, 
orthogonal states | tp^) (m = 2,3, •••,L) together with a m (m = 1,2, •••,L) and 
(3 m (m = l,2,...,L-l), 

a m =(^ m )|Q(0)|^ (m) ) 
/5m = 



-^—{Q (« m _i) | ^ (m - 1} ) - p m -2 1 ^ (m ~ 2) )} , 

Pm-1 



(4) 

(5) 
(6) 



where we define (3$ = and Q (a) = H — a with the Hamiltonian if, so that we obtain 
the tridiagonal matrix 



>1l = 





A 





A 




fa 







a 3 






















































Pl-i 
a L 



(7) 



Pl-i 

The L-th approximate eigenvalue for the ground state is given by the lowest eigenvalue 
of Al- Let us denote the lowest eigenvalue of Al by oll and its eigenvector by 



u 



(L) 



U 



1 ) u 2 ' 



U 



(L) 



l T 



hereafter. Once we know Al, «l and for some 
L, we can in principle evaluate (5l and cvl+i from relations 

^ L) \Q{a L f\^) = {u^YPl, 
(j,™ | Q (a L ) Q (0) Q (a L ) \ j™) = {t#>} 2 /?£ aL+1 , 

with a state 



(8) 
(9) 



u 



(10) 



m=l 
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This would lead us to a larger matrix Al+i and its lowest eigenvalue would give 
us a better estimate of the ground state energy. 

In order to perform a successful numerical evaluation with a small value of L, 
it is necessary to make | if)"') as good as possible. Remember that it is difficult to 
directly calculate a m and (3 m for large systems with such a | ^ }> because our available 
computer memory resources would not be enough to keep whole elements of H \ ip^). 
We therefore truncate each state by operating an M{ v } to it. Details of our stochastic 
selection will be mentioned in section HJ 

3. Trial state 

Before describing the stochastic selection in our calculation, it would be necessary to 
explain how we prepare the trial state | tp^) to start with. In this section and in the 
next section we concentrate our attention on the N s = 48 case, which involves more 
technical issues than the N s = 27 case does. 

First let us comment on our basis by which we describe the states | ip( m '). The 
representation we use is the conventional one where a state is represented by z- 
components of all spins of the system. Here we add an assumption on symmetries. On 
a triangular lattice with 48 spins there are 48 translation operators as well as 6 rotation 
and 2 inversion ones that commute with the Hamiltonian Since we expect that the 
ground state is invariant under these operations [IB], we construct a basis {| <E>j)} which 
ensures all of these invariances. Note that each basis state | <3>j) therefore contains 
maximally 576 (= 48 x 6 x 2) degenerate spin configurations in it. Total number of 
| $j) with S z = therefore is about 6 x 10 10 . 

Now we come to the starting point of our numerical work, which is to calculate 
coefficients of | r/?W) with the basis stated above. We do it in the same way that we 
did in our previous work [T^], where we obtained an approximate ground state | tp\) for 
the spin system on a 36-site triangular lattice. The only difference is that we include as 
many degenerate Ising-like configurations as possible in the initial trial state this time. 
The basic idea for this improvement comes from the Wannier's rigorous proof[24 that 
a classical antiferromagnetic Ising system on a triangular lattice is heavily degenerated 
in its minimum energy, namely its energy at zero temperature, which is —N s /A for the 
iVs-site system. 

Using the conventional Monte Carlo method at low temperature (T = 0.5), where 
the classical energy is used as the Boltzmann weight, we pick up states with the classical 
minimum energy, which is —12 for the N s =48-site system. We find 13 087 | $j)'s 
that fulfill conditions ($j | H | $j) = —12.0 and S z = 0. We then form an initial 
trial state | \l/o) by linearly combining all of them with the equal weight 1/V43087. 
Within this partial Hilbert space with 13 087 basis states, we next pursue the state | \l/ t ) 
which has the lowest energy in the conventional exact diagonalization. We observe that 
(* t | H | %) = -20.1. 

The final stage to calculate a trial state | ^ 1 >) is to repeat following procedures 
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until the obtained value does not change within five decimal digits. 

(i) Extend the partial Hilbert space {| $i) ; ($j | H | $;) = -12.0, S z = 0} by 
operating H a few times to | \&t) until the available computer memory is exhausted. 
The maximum number of basis states we can permit is ~ 1 x 10 8 . 

(ii) Within the Hilbert space determined in (i), pursue the state | with which 

| H | ty' t ) is as low as possible. 

(iii) Form a state | \&") by dropping small coefficients of | We usually request that 
the size of the reduced Hilbert space, which is spanned by basis states included in 
| should be a few percents of the one obtained in (i). 

(iv) Replace | * t ) by | 

After this calculation we obtain the trial state | ip^ 1 ') from the last | We observe 
that | comprises 75 746 657 basis states. Then we obtain ot\ = (ip^ \ H \ ip^) = 
—25.950. It should be noted that we can calculate this inner product exactly because 
we do not need to keep the 'outer' part of H \ ip^) which are orthogonal to the vector 
space attached to | ipW). 

4. Calculations 

With the Hamiltonian and the trial state | ip^) described in the previous section, we 
calculate the matrix elements in Al up to L = 4. This section is to show in detail how 
we carry out the calculations with stochastic selections. Note that all the values and 
the conditions stated below in this section are those for the N s = 48 lattice. Comments 
on the results summarized in tables Q and El will be presented in the next section. 

First step here is to estimate /?i and a 2 so that we can solve the eigenvalue 
problem with A%. Using the recursive SSS |T3j . we calculate the most probable value of 
(3\ = (iffO-) | Q ( ai f | from the statistical average of 

<</> (1) | ■ Q (at) M H 2 )} Q (a,) M {r)W] \ . (11) 

Here we insert a symbol • after (ip^ \ in order to represent that we calculate the inner 
product between | ip^) and Q (cti) M^^yQ («i) M^wy \ ip^). Note that each random 
choice operator in the recursive SSS method depends on the preceding intermediate 
state. We generate each M{ v } adjusting the value of e to be as small as possible 
for a Pentium IV machine equipped with a 2 Giga byte memory. Our result from 
22 samples is (3\ = 1.7082 ± 0.0021, where the error is estimated by the standard 
deviation of the data. When we evaluate a 2 , we approximate | ip^) = Q (cti) | ip^ 1 '}/^ 
by Q M{ v } | ip Namely we calculate 

| ■ Q (a x ) M m Q (0) M {v(2)} Q (a,) M {rjW} | (12) 

to obtain (3\a2- The statistical average from 40 samples gives us a 2 — —12.007 ± 0.066. 
Here we take account of both the statistical error for (fT2"|) and the error from (3\. It is 
easy to find that a 2 = -26.1559, uf ] = 0.99279 and = -0.11986 with the above 
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values of a%, Pi and a 2 . Error estimations on these quantities are carried out as follows. 

~ (2) (2) 

Let / be a 2 , u\ or u 2 . We assume, with errors APi and Aa 2 , 



/(ai, Pi ± Ap\, a 2 ± Aa 2 ) ~ /(ai, p\, a 2 ) ± 

9/ 



and numerically calculate 
matrices 

A ± Aft 



ft ± A/3i 

Q' 2 



and 



df_ 

dPx 

Of 



AA 



+ 



9/ 
<9a 2 



Aa 2 



(13) 



da- 



Aa 2 



by solving eigen problems of 



«i Pi 
Pi a 2 ± Aa 2 



Then we proceed to estimate p 2 and a 3 using (jSJ) and Q with L 



approximate | ip^) 



(2) 



tk) + u { 2 ] | i/> 3 ) by P 2 (M w ) | V (1) >, where 



ft M. 



Pi 



(2) 



-Q(ai)M iv} . 



We calculate statistical averages of 

| • P 2 (M {v m } ) Q(a 2 )M H6)} Q(a 2 )M {vi s )} P 2 (m {v w } ) | ^ (1) > 



(14) 
2. We 

(15) 
(16) 



and 



| ■ A feijW}) Q («2) M {J7 (8 )} Q (0) M {r , (6)} Q (a 2 ) M { „( 5)} P 2 (m { „w>) | (17) 



to evaluate {u^} 2 P 2 and {ii 2 2 ^ } 2 /3f 0:3, respectively. From 34 and 128 samples of these 
quantities, we estimate p 2 and a$. We also evaluate their errors, taking the error from 



u 2 into account. Calculations for a 3 and Urn S 3X6 then straightforward including error 
estimations. For the results, see table El 

In estimations of ft and «4 we calculate inner products between 



(3), 



Q (a 2 ) M {ri(12)} P 3 (M { 



If,®) 



and 



R [M {r)( iti )} j Q (a 2 ) M {ri (i 5 ) } P 3 ^M { ^(i4) } , M {vi i S)} 
where we use abbreviations 



0«) 



PzlM {v} ,M W} 



u 



(3) Pi (3) 



,(3) 



U 2 u 3 



(3) 



, + ^Q(a 2 )M w }>Q(e* 1 )M {?/} 



R[M. 



1 to calculate {u^} 2 Pl 

Q (0) M{^} to calculate {u^YPlon 



(19) 



(20) 



(21) 



for convenience, so that we can generate as many samples as possible within a limited 
CPU time. The results for ft and 0:4, which are obtained from 6240 and 4069 samples 
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respectively, are in table 121 The maximum number of the selected basis states amounts 
to ~ 1.1 x 10 8 . J 

Finally, one comment would be necessary from technical point of view. It should 
be noted that there are variety of ways to calculate samples using the recursive SSS 
method for the following two reasons. 

(i) Suppose we calculate 0f, for example. For this purpose we calculated inner products 
(JTTJ). Yet, the statistical average of inner products between Q (a\) Mum | ip^) and 
Q («i) Mj-^17)} | ip^) will also give us the same quantity. 

(ii) Locations of random choice operators are not uniquely determined. For instance, 
it is possible to employ 



(2) 

U2 QK 



Pi 



M. 



and 



pUm {v} ,m W} 



u 



(3) 



- ^u {3) + 



u 



(3) 



01 



+ 



(•3) 



Q (a 2 ) > M {n} Q (at) 



instead of P 2 ( M^j ) and P-. 



M {ri} ,M W} 



0102 

respectively 



(22) 



M W} , (23) 



Although statistical averages obtained in these ways will theoretically agree with each 
others, their standard deviations might be different. It is not a priori clear what way is 
best for numerical calculations. In order to present an example which shows how much 
difference is actually observed, we calculate 300 samples of 



| .i* (M {}7 (23) } , M { „(22) } J Q (a 3 ) M W 2 1)} Q (a 3 ) M {r)(3 o )} P^ [M^y, M H is )} j | , (24) 
whose statistical average also should give 

{u i 3 ) Y0i- The result is 3 = 10.8 ±11.0 
with a fixed value of e = 0.0025, and we observe almost the same variance as in the 
measurement with P 3 (M^y, M^>y 



5. Summary and discussions 

In this paper we combine our recursive SSS method with the Lanczos approach so 
that we can estimate the ground state energy of the spin-1/2 quantum Heisenberg 
antiferromagnet on a 48-site triangular lattice. 

In order to examine whether the method works well, we study the 27-site system 
for which the ground state energy is exactly known to be E/JN^ = —0.1867404. Our 
results for the 27-site system in tableHJgive a satisfying upper limit E / JN^ < a^/ JNb = 
—0.18506 ± 0.00037. It should be noted that we obtain this upper limit within a limited 
partial Hilbert space. 

| Similar calculations for @4 and would be possible with a swifter computer. We estimate that the 
CPU time necessary to calculate these quantities with our Pentium IV machine is about 300 times as 
long as the CPU time we spent for (3% and a,\. 
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Now we summarize our results for the 48-site system. Starting from a state | ip^) 
with ~ 7.6 x 10 7 basis states, each of which being a representative of maximally 576 
translation-, rotation- and inversion- invariant configurations, we successfully calculate 
elements of the tridiagonal matrix A4, namely a\, a 2 , a 3 , a 4 , /3i, j3 2 and j3 3 presented 
in table El Our best estimate for the upper bound of the ground state energy is given 
by a 4 . The result per bond is <S 4 / JN b = —0.1833 ± 0.0003. 

This value should be compared with values obtained by other methods. In a figure 
presented by Capriotti et. al. who study the system using the Green function 
Monte Carlo method, we see that E/JN b ~ —0.185 when N s = 48. The variational 
Monte Carlo study [23] presents a value Ej JN b = —0.185 ± 0.001 for the 48-site system. 
Richter et. al. [4 , on the other hand, made the finite-size extrapolation using the 
results for N s = 24, 30 and 36 obtained by the exact diagonalization method. Using 
the scaling formula eo (v^s) — A + A 3 / (v^Vs) 3 + ® (^T 2 ) f° r the lowest energy per 
site e = E/JN S = 3E/JN b , they obtain A /3 = e (oo)/3 = -0.1842.§ The value 
for the 48-site system in this formula, which we calculate from this Aq and the value 
e (v^6)/3 = -0.1867912 [16;, is e {V48)/3 = -0.1859. Thus we see that our result 
is consistent with those obtained by the Green function Monte Carlo method and the 
variational Monte Carlo method. This indicates that the recursive SSS method combined 
with the Lanczos approach proves to be one of techniques applicable to frustrated 
quantum systems. 

A few remarks are in order. 

Applications of the SSS method in study of other models, especially of frustrated 
ones, are in prospect. For example, we have studied the 64-site Shastry-Sutherland 
model combining the simple SSS method with the power method [T2]. where we obtained 
results to indicate that the intermediate phase exists. Much improved results on this 
model are expected with our new method presented in this paper. 

The result a^/JN b = —0.1833 ± 0.0003 on the 48-site triangular lattice is much 

more precise compared with our previous result E^ t /JN b = —0.1856 ± " on the 

36-site lattice jT3] - One reason for this improvement is that, instead of the power method, 
we adopt a Lanczos approach with which the fitting procedure is not necessary. Another 
reason is that we improve the trial state noticing the existence of heavily degenerate 
Ising-like states that have the lowest energy in the classical Ising system on the triangular 
lattice. This idea, inspired by the Wannier's proof[2I] on the classical Ising system, is 
encouraged by observations in numerical study of quantum systems on small lattices. 
For N s = 12, 21 and 27 lattices we observe that a large part of the ground state is 
formed by those degenerate Ising-like basis states. 

The final remark is on the merits of the SSS method. First of all, we emphasize 
that this method is mathematically justified. It is guaranteed that expectation values 
for any operator and any state are given correctly by the statistical averages. Secondly 

§ Simple extrapolation using our results on N s = 36 and 48 lattices yields an upper bound 
E/JN b < -0.1790 ± 0.0025 in the N s -> 00 limit. 
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our method is quite general in a sense that the SSS method can be substituted for the 
importance samplings in usual Monte Carlo methods 12 • This should be compared 
with the variational techniques, for which one needs deep insights on properties of the 
system under study[20] |2T]. Thirdly one can easily join the SSS technique to various 
methods such as the power method, the Lanczos method and the variational method. 
We expect this technique will be helpful to improve many methods used in numerical 
studies. 
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Table 1. Estimated values of as and /3s in the Lanczos approach for the N s = 27 
system, to which we know that E cxact = — 15.125972^^. For this lattice we simply 
collect all S z — 1/2 states to form a basis. The size of the full Hilbert space is 
therefore 27Ci3=20058 300. We start from a trial state with 1072 935 basis states 
obtained after the procedures in section |3 for which ot\ = —14.5680. The last three 
rows in the table present exact values which are calculated in the full Hilbert space 
without any stochastic selections. When we calculate as and /3s inside a limited partial 
Hilbert space, we set the maximum size of the space to be 4 500000. 



L 


2 


3 


4 




1.5405 ± 0.0015 


6.14 ± 0.18 


6.09 ± 0.55 


ClL 


-4.999 ± 0.014 


-3.53 ± 0.39 


-1.8 ± 1.5 


a L 


-14.8099 ± 0.0023 


-14.927 ± 0.010 


-14.990 ± 0.030 


Pl-i 


1.54010 


6.23049 


6.06325 


^exact 


-5.01580 


-3.72843 


-1.85916 


^exact 


-14.8101 


-14.9354 


-15.0061 



Table 2. Estimated values of as and /3s in the Lanczos approach for the N s = 48 
system. The maximum number of basis states is limited to 1.1 x 10 8 . As is described 
in section|21we first calculate a trial state with 75 746 657 basis states, for which a± is 
—25.950. Then, following equations (|llf> and (|12|) in section^] we estimate /3i and a.2- 

~ (2) (2) 

Thirdly we calculate aa and u\ by an exact diagonalization. In a similar way to 

that stated in section 0] in detail, we evaluate quantities for L = 3 and L = 4. 



L 


2 


3 


4 


(3l-i 


1.7082 ± 0.0021 


12.63 ± 0.16 


13.2 ± 2.0 


a L 


-12.007 ± 0.066 


-1.43 ± 0.45 


12.1 ± 6.6 


a L 


-26.1559 ± 0.0011 


-26.3189 ± 0.0096 


-26.393 ± 0.049 


(L) 
u\ 


+0.99279 ± 0.00007 


+0.9719 ± 0.0016 


+0.955 ± 0.013 


(L) 
u\ 


-0.11986 ± 0.00056 


-0.2101 ± 0.0051 


-0.248 ± 0.024 


(L) 




+0.1066 ± 0.0048 


+0.153 ± 0.030 








-0.052 ± 0.024 



